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INTRODUCTION 

This report discusses the progress made under Crustal Dynamics Project funding by the 
University of Texas Center for Space Research on a wide variety of topics, including geodesy, 
geodynamics, and satellite dynamics. The report is derived from a manuscript submitted to the 
American Geopohysical Union Monograph dedicated to the Crustal Dynamics Project. 

The ability to range accurately from the Earth’s surface to satellites carrying laser cube-comer 
reflectors, along with the launching of the Lageos satellite in May 1976, has provided a unique 
capability for studying global solid Earth dynamics [Johnson et al., 1976]. During the period 
between Lageos launch through completion of the Crustal Dynamics Project at the end of 1991, the 
satellite laser ranging (SLR) technique has evolved into one of the fundamental geophysical and 
geodetic measurement techniques [Tapley et al., 1985; 1990, Frey and Bosworth, 1988], The 
primary goals of the SLR system development and demonstration were to measure tectonic motion 
and Earth orientation. In the early 1980’s, the technique achieved operational status in measuring 
Earth orientation and global baselines with the precision required to study plate motion and 
deformation. In addition, the technique has demonstrated unique abilities to measure both the 
constant and time-varying gravitational field properties of the Earth, to provide a unique terrestrial 
reference frame tied to the Earth’s center of mass, and to study the dynamical effects of general 
relativity. It is currently regarded as a required tracking system for altimetric satellites such as 
TOPEX/POSEIDON and ERS-1. 

The Lageos satellite, spherical in shape and covered with 426 comer cube reflectors, 422 of fused 
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silica for reflectance in the visible spectrum and 4 of germanium for infrared reflectance, was placed 
in a nearly circular orbit with an altitude of one Earth radius [Johnson et al., 1976; Cohen and Smith, 
1985]. The combination of satellite design and orbit has provided an ideal ranging target for the 
SLR systems, which were originally developed in the 1960’s, and have developed from meter level 
precision to the current subcentimeter level. The Lageos orbit provides a stable reference frame for 
studying the rotation of the Earth, the relative motion of points on the Earth’s crust, and the time 
variation of the long- wavelength component of the Earth’s gravity field. The Lageos orbit has a 
mean semimajor axis of 12,271 km, and eccentricity of 0.0044, and an inclination of 109.84° with 
respect to the Earth’s equator. The orbit plane completes one rotation with respect to the true-of-date 
equinox in 1050 days in a prograde sense, while the perigee completes one revolution with respect to 
the Earth’s equator in 1680 days in a retrograde sense. The satellite’s altitude reduces the effects of 
uncertainties in the models for the high degree and order portion of the Earth’s gravitational field and 
the effects of atmospheric drag. In addition, the satellite has a beryllium-copper core to increase its 
mass (407 kg) and a 60 cm diameter, thus giving the satellite a very small area-to-mass ratio of 
6.95x10 m /kg. This further reduces the effect of difficult-to-model surface forces such as neutral 
and charged particle drag, radiation pressure, and thermal forces. The resulting orbit stability give 
the Lageos satellite a projected lifetime of over 500,000 years. Significant periods associated with 
the Lageos orbit are provided in Table 1. 


Table 1. Significant Lageos Orbit Characteristics 

Parameter 

Period (days) 

Orbital period 

Node period wrt inertial space 
Node period wrt Sun 
Perigee period wrt inertial space 

0.1566 (225 minutes) 
1050 (prograde) 

560 (prograde) 

1 680 (retrograde) 
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THE SLR GROUND NETWORK AND DATA SET 


In the investigation described here, the University of Texas orbit analysis system, UTOPIA, was 
used to obtain a dynamically consistent solution to 15.1 years of Lageos laser range data. The data 
set spanned the period from May 7, 1976 to July 3, 1991, and included 561,708 normal points, 
constructed from two- and three-minute averages of the laser range observations from 117 different 
SLR sites. The single-shot ranging precision of the various systems varied from about a meter for 
some systems in 1976 to better than 10 mm for many of the systems at the end of the solution 
interval. The normal point precision is about an order of magnitude better. An important component 
of the tracking system coordinates is the offset, or eccentricity, between the laser site and the 
appropriate geodetic monument. These eccentricities are measured using ground surveying 
techniques with accuracies better than 1 cm. The eccentricity file used by UTCSR for this study was 
compared against a similar file maintained at the Goddard Spaceflight Center and indicated no 
significant discrepancies [Sellars, 1989]. The nominal epoch station coordinates used for this 
investigation were Lageos derived coordinates referred to as SSC(CSR)91L02. The solution 
approach is described by Watkins [1990], The nominal plate motion model was the no net rotation 
absolute motion model AMO-2 of Minster and Jordan [1978]. The epoch for the plate motion was 
1988.0. 

As shown in Table 1, Lageos has a period of 225 minutes, providing 3-5 passes per day for most 
locations on the Earth. However, most stations attempt to acquire ranges for only the passes in a 
single 8-hour shift, which is usually scheduled during local nighttime. To reduce the computational 
burden, average or normal points are formed from the full-rate data collected by the station. These f 
normal points are essentially average ranges over a selected duration, such as 2 or 3 minutes, and 
have been demonstrated to retain the geophysically and geodetically useful information of the full- 
rate data set [Smith et ah, 1985; Tapley et ah, 1985]. Normal points are used not only to relieve 


1.2 



computational effort but also to reduce the random or white noise content of the laser range 
measurements. The procedure accomplishes this by averaging a large number of raw ranges, and 
thus will not eliminate systematic errors which may be relatively constant or slowly varying over the 
two- or three-minute normal point window. The formation of the normal points used for this study 
follow the Herstmonceux recommendation outlined at the Fifth International Workshop (1984) on 
Laser Ranging Instrumentation. As Table 2 indicates, this approach has led to the set of 561,708 
normal points which were used in this study. The mean precision for this set of 1 17 sites was 7.3 cm 
when averaged over the entire span, although many of the recent sites have average precisions of less 
than 1 cm. 

Data Weighting 

The quality and quantity of the Lageos range data has varied widely during the 15 years of the 
mission. The variation in data quality requires that a complex weighting algorithm be applied if the 
best estimates of the satellite orbit and geodetic parameters are to be obtained. 

The data quality variation can be thought of as being of two types. The first type is the gradual 
increase in data quality with time as the laser instrumentation has improved. This is demonstrated 
by Figure 1, which plots the rms fit of normal points in 15-day arcs during the first 10 years of the 
mission, using the models described in the next sections. The resulting curve in Figure 1 was 
approximated with a four-part linear model with the following node points: 


1976 (MJD 42905): 

70 cm 

1979 (MJD 44162): 

40 cm 

1983 (MJD 45578): 

20 cm 

1986 (MJD 47000): 

12 cm 


The weight of each normal point is linearly interpolated between the node points and is constant at 
the value of the last node point for times greater than the time of the last node. 
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The second type of variation in data quality is the variation from station to station within the 


network due to the variety of hardware systems in use at the tracking stations. This inhomogeneity 
is demonstrated by the results in Table 2. This table provides an estimate of the internal precision of 
each station in the network (column titled Prec. Est.), averaged over all data retained in the solution. 
This variation is modeled as a station dependent noise level which is added in an rms sense to the 
linear weighting described above. This correction was used only for those stations whose noise level 
is substantially higher than the network average at the time of operation of the site. In addition, 
recently occupied mobile sites are assigned slightly higher noise levels to reflect errors in nominal 
site positions. For this study, the systems were assigned the incremental noise levels in centimeters 
indicated in Table 2 (column titled Inc. Sig.). Stations with 0 in this column receive only the 
weighting assigned through the time dependent algorithm described above. 


Table 2. SLR Station Performance 
May 1976- August 1991 



Station 

No. of 
Passes 

No. of 
Obs. 

Prec. Est 
(cm) 

Inc. Sig. 
cm 

1181 

Pottsdam, Ger. 

559 

4048 

9.1 

30 

1873 

Simeiz, Uk. 

82 

590 

7.3 

100 

1884 

Riga, Lat. 

121 

794 

6.9 

100 

1893 

Katsively, Uk 

41 

304 

4.6 

100 

1953 

Santiago de Cuba, Cuba 

106 

532 

8.4 

100 

7035 

Otay Mt., USA 

48 

717 

0.5 

0 

7046 

Bear Lake, USA 

45 

651 

3.6 

10 

7051 

Quincy, USA 

137 

1271 

6.0 

0 

7062 

Otay Mt., USA 

265 

1958 

2.9 

0 

7063 

STALAS, GSFC, USA 

446 

4219 

4.2 

0 

7065 

GSFC, USA 

3 

19 

8.3 

0 

7067 

Bermuda 

29 

161 

3.1 

0 

7068 

Grand Turk, Bahamas 

4 

22 

8.4 

0 

7080 

McDonald Obs., USA 

929 

11369 

1.0 

0 

7082 

Bear Lake, USA 

117 

843 

4.1 

0 

7084 

Owens Val., USA 

22 

152 

14.3 

0 

7085 

Goldstone, USA 

20 

135 

8.6 

0 

7086 

Mcdonald Obs., USA 

1270 

13197 

2.3 

0 

7090 

Yaragadee, Aust. 

3661 

45811 

1.9 

0 

7091 

Haystack Obs., USA 

412 

3860 

5.5 

0 

7092 

Kwajalein 

55 

497 

9.2 

0 

7096 

American Samoa 

124 

953 

6.7 

0 

7097 

Easter Is., Ch. 

198 

2398 

1.2 

10 

7100 

GSFC, USA 

5 

41 

6.7 

0 

7101 

GSFC, USA 

9 

67 

7.2 

0 
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7102 

GSFC, USA 

229 

2019 

6.7 

0 

7103 

GSFC, USA 

60 

578 

4.8 

0 

7104 

GSFC, USA 

13 

110 

3.1 

0 

7105 

GSFC, USA 

1995 

24343 

0.8 

0 

7109 

Quincy, USA 

2652 

38505 

1.2 

0 

7110 

Mon. Peak, USA 

3129 

41504 

1.2 

0 

7112 

Platteville, USA 

624 

6427 

3.5 

0 

7114 

Owens VaL, USA 

373 

3316 

3.5 

0 

7115 

Goldstone, USA 

379 

3293 

2.7 

0 

7120 

Maui, USA 

368 

3460 

2.1 

0 

7121 

Huahine, Fr. Poly. 

262 

2125 

3.1 

0 

7122 

Mazatlan, Mex. 

1356 

17573 

1.4 

0 

7123 

Huahine, Fr. Poly. 

179 

2232 

0.9 

0 

7125 

GSFC, USA 

35 

325 

1.7 

0 

7130 

GSFC, USA 

30 

253 

1.0 

0 

7210 

Maui, USA 

2084 

23249 

1.1 

0 

7220 

Mon. Peak, USA 

63 

511 

4.0 

0 

7236 

Wuhan, PRC 

52 

678 

2.6 

60 

7265 

Mojave, USA 

48 

438 

1.9 

0 

7274 

Mon. Peak, USA 

48 

384 

3.9 

0 

7288 

Mojave, USA 

185 

2822 

0.6 

10 

7295 

Richmond, USA 

86 

1029 

1.6 

10 

7400 

Santiago, Ch. 

42 

339 

2.3 

10 

7401 

Cerra Tolola 

222 

3008 

1.2 

10 

7403 

Arequipa, Peru 

150 

2174 

1.0 

10 

7510 

Askites, Gr. 

258 

2589 

2.0 

10 

7512 

Rhodes, Gr. 

176 

1921 

2.0 

10 

7515 

Dionysos, Gr. 

247 

2969 

2.1 

10 

7517 

Roumeli, Gr. 

225 

2604 

1.3 

10 

7520 

Karitsa, Gr. 

95 

1052 

2.1 

10 

7525 

Xrisokellaria, Gr. 

126 

1570 

0.9 

10 

7530 

Bar Giyyora, Is. 

97 

920 

4.5 

30 

7540 

Matera, It. 

32 

225 

1.4 

10 

7541 

Matera, It. 

61 

535 

1.9 

10 

7543 

Noto, It. 

49 

529 

1.4 

10 

7544 

Lampedusa, It. 

145 

1577 

2.0 

10 

7545 

Cagliari, It. 

120 

1323 

1.6 

10 

7546 

Medicina, It. 

11 

85 

2.3 

10 

7550 

Bassovizza, It. 

49 

257 

2.5 

10 

7575 

Diyarbakir, Tur. 

89 

1183 

1.6 

10 

7580 

Melengiclic, Tur. 

91 

1227 

0.4 

10 

7585 

Yozgat, Tur. 

84 

909 

2.0 

10 

7587 

Yigilca, Tur. 

83 

1102 

0.6 

10 

7596 

Wettzell, Ger. 

37 

262 

2.1 

0 

7597 

Wettzell, Ger. 

10 

75 

1.0 

0 

7599 

Wettzell, Ger. 

12 

77 

2.1 

0 

7602 

Tromso, Nor. 

45 

405 

2.7 

10 

7805 

Metsahovi, Fin. 

150 

320 

14.1 

100 

7810 

Zimmerwald, Switz. 

828 

10510 

2.3 

10 

7811 

Borowicz, Pol. 

99 

639 

7.3 

100 

7831 

Helwan, Egy. 

198 

1794 

1.7 

30 

7833 

Kootwijk, Neth. 

338 

2399 

11.7 

30 

7834 

Wettzell, Ger. 

1470 

12809 

1.8 

10 

7835 

Grasse, Fr. 

2717 

41961 

2.4 

0 

7837 

Shanghai, PRC 

256 

2251 

5.4 

60 

7838 

Simosato, Jap. 

1343 

13857 

2.9 

10 
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7839 

Graz, Aus. 

1098 

12992 

0.9 

0 

7840 

Greenwich Obs., UK. 

3364 

37127 

1.8 

0 

7843 

Orroral Val., Aust. 

1535 

15321 

2.0 

5 

7853 

Owens Val., USA 

105 

1743 

2.0 

0 

7882 

Cabo San Lucas, Mex. 

54 

641 

0.6 

0 

7883 

Ensenada, Mex. 

35 

273 

0.6 

0 

7885 

Mcdonald Obs., USA 

34 

218 

1.6 

0 

7886 

Quincy, USA 

107 

989 

1.7 

0 

7888 

Mt. Hopkins, USA 

30 

231 

2.5 

0 

7891 

Flagstaff, USA 

36 

281 

2.8 

0 

7892 

Vernal, USA 

66 

507 

4.9 

0 

7894 

Yuma, USA 

45 

221 

2.6 

0 

7896 

Pasadena, USA 

66 

516 

3.1 

0 

7899 

GSFC, USA 

28 

177 

4.2 

0 

7907 

Arequipa, USA 

4046 

46849 

14.9 

30 

7918 

GSFC, USA 

53 

619 

0.4 

100 

7919 

GSFC, USA 

6 

42 

4.9 

100 

7920 

GSFC, USA 

25 

230 

0.5 

100 

7921 

Mt. Hopkins, USA 

686 

6739 

33.0 

50 

7929 

Natal, Braz. 

353 

2152 

31.1 

50 

7939 

Matera, It. 

2242 

32328 

5.3 

20 

7940 

Dionysos, Gr. 

3 

15 

9.7 

100 

7943 

Orroral Val., Aust. 

1381 

14756 

21.2 

50 

8833 

Kootwijk, Neth. 

106 

1148 

2.0 

10 

8834 

Wettzell, Ger, 

39 

349 

1.4 

10 


TOTALS 

48456 

561708 
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METHOD OF SOLUTION 


Numerical Integration 

The solution for the Lageos ephemeris, dynamic model parameters, and geodetic parameters was 
obtained by a weighted least squares batch estimation procedure which requires the solution of the 
differential equations governing the satellite's motion, as well as the state transition matrix for a 
defined period of satellite motion [Tapley, 1973]. For short arc solutions, this can be accomplished 
to a reasonable accuracy by numerically integrating the second order differential equations of the 
orbital motion in the form (Cowell's method) 


? = -A? +?«,?,?) CD 

r 

where ~f represents the perturbing forces to the two body motion. For the longer arcs described in 
this investigation, a modified Encke method was required to numerically integrate Eq. (1) [Lundberg 
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et al., 1991]. The numerical integration method for both the standard Cowell and modified Encke 
approaches was a Krogh-Shampine-Gordon 14 th order, fixed- step integrator described in Lundberg 
[1985]. The step size used was 300 seconds. 

In UTOPIA, the dynamical equations which govern the motion of Lageos were expressed in an 
Earth-centered (non-rotating) cartesian coordinate frame, defined by the mean equator and equinox 
of epoch 2000.0 (J2000.0). This system is realized through the use of the JPL DE-200 planetary 
ephemeris, the 1976 International Astronomical Union (IAU) precession and the 1980 IAU Wahr 
nutation model. Corrections to the VLBI-determined IAU precession and nutation models [Herring, 
1988] were also adopted. The orientation of the tracking stations (the body-fixed frame) relative to 
the celestial ephemeris pole of the Earth are provided through the use of the EOP(CSR)91L02 Earth 
orientation series [Eanes et al., 1991], 

Long- and Short-Arc Solution Procedure 

The determination of geophysical and geodetic parameters using Lageos laser range 
measurements at UTCSR involves a combination of long- and short- arc techniques [Tapley et al., 
1985], The long arc provides the starting point and is used primarily to study long-period 
perturbations in the Lageos orbit. The short arcs are constructed from the residuals of the long-arc 
orbit and form the basis for most geodetic work. A description of both techniques is provided in the 
specific context of the solutions determined for this study. 

Long-Arc Solution 

Using the data set described in Table 2 and a complete force and measurement model, a single, 
dynamically consistent trajectory was fit over the 15. 1-year period from May 1976 through July 
1991. The force and measurement models adhered largely to the IERS Standards, with the exception 
of the use of the TEG-2 mean gravity field developed at UTCSR, and the modelling of a signifcantly 
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more complete ocean and solid tide model, including ocean tide perturbations from over 80 
constituents expanded to degree and order 20 where necessary, and third order solid tides [Tapley et 
al., 1991; Eanes et al., 1991; Casotto, 1989], To achieve an accurate solution over this extended 
time interval, stringent demands must be met on the accuracy of the force models and computational 
software. The numerical integration process involves over 1.6 million steps in over 35,000 orbital 
revolutions, equivalent computationally to a 2700 year integration of the orbit of the Moon, a much 
more conservative dynamical system. 

After converging the orbit through the entire data span in this manner, adjusting only the single 
set of Lageos initial conditions and 15 day along track accelerations, the range residual rms was 
1.28 meters. The residuals from the long-arc solution were mapped into orbit elements using the 
UTCSR software package ELPSOL. A spectral analysis was carried out on the residual orbit 
element time series, and candidates for the sources of errors were identified and added to the 
adjusted parameter list on subsequent iterations to reduce the range residual rms. 


Table 3. Contributions to 1.28 Meter Range 
Residual RMS By Individual Orbit Element 

Element 

Residual 
rms (m) 

Semimajor Axis 

0.002 

Eccentricity 

0.522 

Inclination 

0.127 

Node 

1.211 

Vi (co + m) 

0.498 

Vi (co - m) 

0.477 


The dominant sources of error for each element shown in Table 3, were : semimajor axis - along- 
track acceleration variations with periods shorter than 15 days; eccentricity - odd zonal harmonics 
(constant + variability due to ocean tides and meteorological effects), and thermal surface effects 
(solar Yarkovsky + asymmetric albedo of spacecraft); inclination - ocean tide error dominated by K \ 
and S 2 \ node - even zonal harmonics, particularly the secular variation j 2 , and periodic variability at 
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18.6 years and annual periods; periapse - odd zonal harmonic error (constant + variability), and 
thermal effects; along-track - variations in along-track "drag" as in semimajor axis. The linear 
combinations of periapse longitude and mean anomaly are used because they are more well suited 
for near circular orbits such as that of Lageos. 

The second step in the long-arc procedure uses the 1.28 m orbit obtained by estimating only a 
small subset of parameters to tune a set of force model parameters in order to reduce the range 
residual rms. After adjusting the parameters described in Table 4 (for estimated ocean tides see 
Table 6), the range residual rms was reduced to 28 cm. It should be noted that solar reflectivity was 
not adjusted as a sub-arc parameter in the manner of the along-track acceleration, since previous 
studies at the Center for Space Research have indicated that the accuracy of the estimate of the 
reflectivity does not reach the 0.1% level of the hypothesized solar constant variability unless it is 
adjusted in intervals spanning several years [Willson and Hudson, 1988; Tapley et al., 1989]. This is 
due to the need for the shadowing function to decorrelate the powerful reflectivity parameter from 
other adjusted parameters. 


Table 4. Estimated Parameter Summary 

Parameter 

Frequency 

Initial conditions 

Once 

/ 2 >^3 

Once 

h 

Once 

Ocean tides 

Once 

Solar reflectivity 

Once 

Along-track acc. 

15 days 


Short-Arc Solutions 

The 28 cm residuals from the final long arc fit were dominated by small model errors in ocean 
tide perturbations, including seasonal variability and thermal (Yarkovsky) forces [Eanes and 
Watkins, 1991]. Discussions of tide model errors including errors due to omission, commission, and 
to seasonal variability in the Lageos orbit are given by Eanes and Watkins [1991] and Casotto 
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Table 5. Contributions to 0.28 Meter Range 
Residual rms by Individual Orbit Element 

Element 

Residual 
rms (m) 

Semimajor Axis 

0.001 

Eccentricity 

0.271 

Inclination 

0.066 

Node 

0.126 

Vi (co + m) 

0.259 

Vi (oo — m) 

0.232 


[1989]. The rms contribution of errors in each orbit element is presented in Table 5. These errors, 
because of both the magnitude and spectrum, can alias into the solutions for the geodetic parameters. 
The effects of long period force model errors must be removed from the residuals before an accurate 
solution for geodetic parameters can be obtained. The short-arc solution is designed to achieve this 
result, since the long period error can be accomodated in the estimate for the initial conditions in a 
manner similar to the classical variation of parameters approach used in celestial mechanics. The 
length of the short arcs were dependent on the data density and varied from 15 to 3 days according 
to: 


1976 (MJD 42905) - 1979 (MJD 44162): 15 day 

1979 (MJD 44162) - 1983 (MJD 45578): 6 day 

1983 (MJD 45578) - 1989 (MJD 47585): 3 day 

These cutoffs for the arc lengths were chosen to make the a posteriori uncertainties on the estimated 

orbit parameters more uniform over the data span [Watkins et ah, 1989]. The Earth orientation 

parameters, station coordinates and other parameters of primary geodetic interest were adjusted 

simultaneously with the short-arc adjustments. 

DYNAMICAL MODEL INVESTIGATIONS 

As mentioned previously, the long arc technique is ideally suited to studying the accuracy of the 
dynamic models used to propagate the satellite orbit. If the dynamic model used to compute the 
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long-arc orbit is imperfect, the tracking data cannot be fit to the measurement precision, and 
systematic range residuals will result. By fitting the residuals with piecewise constant adjustments to 
the classical orbital elements over successive time intervals that are short with respect to the length 
of the long arc, a time series of the error in each orbital element is obtained. Spectral analysis of the 
time series provides a means of identifying what parts of the dynamical model need to be adjusted. 

Gravitational Forces 

Secular and Tidal Variations 

Figure 2 shows an example of the above described process. The top panel shows the Lageos 
inclination residuals from May 1976 to July 1991 as computed using the nominal dynamical model 
with ocean tides fixed to the values from Schwiderski [1980]. Because errors in the dynamical 
model drive the derivatives of the orbital elements, a more direct comparison of the size of the model 
errors is obtained by analyzing the derivative of the orbital element time series. Figure 2c shows the 
derivative of the inclination residuals, the observed inclination excitations, and Figure 2d shows the 
power spectrum of these excitations. Peaks in the power spectrum at periods of 1050 days, 280 days, 
and 14.03 days are marked and labeled K } , S 2 an d ^ 2 « These peaks show that there are errors in the 
nominal ocean tide model for these constituents. In particular this inclination signal points to the 
need to adjust the prograde degree 2 order 1 harmonic of the K x tide and the prograde degree 2 order 
2 harmonics of the 5*2 and M 2 tides. 

Figure 2b shows the inclination residuals after these parameters (and others) have been adjusted 
using the Lageos data. The systematic signals are significantly reduced. The RMS of the inclination 
residuals about the best fitting line is 10 mas in Figure 2a and is 3 mas in Figure 2b. The slope of the 
remaining inclination residuals in Figure 2b is 0.3 mas/yr. The remaining signals, although 
substantially smaller, still indicate that further model improvement is possible. 
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Analysis of the other orbital element residual time series show that many other periodic signals 
exist that can be removed by adjusting tide parameters. Tables 6, 7, and 8 show the results of the 
adjustment of these tide parameters. The phase definition in these tables is adopted from that of 
Schwiderski [1980] as given in the IERS Standards [McCarthy et al., 1989], 

The S a tide, although not included in the 11 constituents computed by Schwiderski, can, 
nevertheless, be estimated from the other long period tides by assuming that the admittance of the 
ocean’s response varies smoothly with frequency [Eanes et al., 1983; Christodoulidis et al., 1985]. 
Since S a is a small term in the tidal potential the resulting estimate for the S a tide parameters are also 
small. The Lageos values, however, are not small and clearly indicate that non-tidal sources of 
dynamical model error are present. 




Table 6. 

Long Period Ocean Tide Solutions 


C20 

e 20 

C 30 

630 

Source 


cm 

deg 

cm 

deg 


2.17 

23 

11.4 

282 

Lageos (this paper) 


0.17 

264 

0.01 

46 

Schwiderski [1980] 


2.83 

39 

1.98 

245 

Starlette, Cheng et al. [1990] 

Ssa 

1.81 

267 

1.98 

81 

Lageos (this paper) 


1.24 

222 

0.06 

2 

Schwiderski [1980] 


1.59 

253 

0.78 

69 

Starlette, Cheng et al. [1990] 


1.29 

262 

0.53 

162 

Lageos (this paper) 


1.06 

259 

0.06 

94 

Schwiderski [1980] 


1.42 

246 

- 

- 

Starlette, Cheng et al. [1990] 

M f 

2.86 

243 

2.61 

281 

Lageos (this paper) 

1.70 

252 

0.19 

148 

Schwiderski [1980] 


2.84 

242 


- 

Starlette, Cheng et al. [1990] 


In the case of C \ 0 , the estimate of the 2-cm S a "ocean tide" is actually the result of seasonal 
redistribution of mass in the atmosphere and hydrosphere [Gutierrez and Wilson, 1987; Cheng et al., 
1989; Chao and Au, 1991], The 2-cm value for C 2 o is equivalent to an annual variation of J 2 with 
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Table 7. 

Diurnal Ocean Tide Solutions 


C Si 

£21 

^31 

e 31 

ex 1 

eJi 

Source 


cm 

deg 

cm 

deg 

cm 

deg 

Q\ 

0.72 

314 

0.08 

297 

0.48 

281 

Lageos (this paper) 


0.54 

314 

0.31 

107 

0.29 

289 

Schwiderski [1980] 


0.72 

295 

- 

- 

0.87 

264 

Starlette, Cheng et al. [1990] 

Ox 

2.46 

308 

0.91 

64 

1.77 

301 

Lageos (this paper) 


2.42 

314 

1.31 

84 

1.43 

276 

Schwiderski [1980] 


2.66 

327 

1.05 

63 

2.25 

296 

Starlette, Cheng et al. [1990] 

p 1 

0.85 

323 

0.47 

48 

- 

- 

Lageos (this paper) 


0.90 

314 

0.30 

40 

0.63 

258 

Schwiderski [1980] 


0.99 

331 

0.86 

1 

0.78 

267 

Starlette, Cheng et al. [1990] 

Si 

0.14 

328 

1.12 

232 

- 

- 

Lageos (this paper) 


0.02 

315 

0.01 

37 

0.01 

256 

Schwiderski [1980] 

K x 

2.47 

326 

1.31 

31 

- 

- 

Lageos (this paper) 


2.81 

315 

0.89 

34 

1.91 

254 

Schwiderski [1980] 


2.68 

325 

1.41 

347 

2.59 

254 

Starlette, Cheng et al. [1990] 


an amplitude of 24 x 10“ n . The S sa tide results also differ substantially from Schwiderski and are 
due to a semiannual variation of the mass distribution in the atmosphere and oceans. On the other 
hand, the 1 1.4 cm amplitude of S a , represented by the value for C 30 , is too large to be explained by 
mass transfer in the atmosphere. Since similar analysis using Starlette [Cheng et al., 1989] does not 
observe the same large signal, this anomaly indicates that there is some nongravitationa! term that 
needs correction in the Lageos dynamical model. In addition to the anomalous value of S a C 30 , the 
value of S\ C 31 obtained using the Lageos data is also large. These two anomalies are probably 
related and are discussed below. 

The Lageos and Starlette results for Mj C \ 0 are both consistently larger than the Schwiderski 
value and show the same phase. The likely explanation is that both estimates are aliased by error in 
the nominal model for the fortnightly variation of UT1 [Yoder et al., 1981]. For Lageos, a 1-cm 
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Table 8. 

Semi-diurnal Ocean Tide Solutions 


C72 

e 22 

C 32 

e 32 

C 42 

£42 

Source 


cm 

deg 

cm 

deg 

cm 

deg 

n 2 

0.72 

325 

0.38 

205 

0.29 

132 

Lageos (this paper) 


0.65 

322 

0.11 

172 

0.21 

142 

Schwiderski [1980] 


0.92 

329 

- 

- 

0.13 

148 

Starlette, Cheng et al. [1990] 

m 2 

3.30 

321.4 

0.27 

155 

1.00 

135 

Lageos (this paper) 


2.96 

310.6 

0.36 

169 

1.00 

125 

Schwiderski [1980] 


3.22 

319.3 

0.12 

161 

1.15 

121 

Starlette, Cheng et al. [1990] 

Si 

0.84 

306 

0.38 

214 

- 

- 

Lageos (this paper) 


0.93 

314 

0.26 

202 

0.37 

103 

Schwiderski [1980] 


0.83 

302 

0.23 

192 

0.32 

89 

Starlette, Cheng et al. [1990] 

K 2 

0.31 

320 

0.19 

244 

- 

- 

Lageos (this paper) 


0.26 

315 

0.09 

195 

0.11 

104 

Schwiderski [1980] 


0.29 

316 

0.08 

243 

0.10 

97 

Starlette, Cheng et al. [1990] 


change in the amplitude of MyCjo corresponds to 0.02 ms or less than 2 percent of the total 
variation of UT1R-UT1. The causes of the large estimates of M m C 20 and My C 30 are not known at 
this time. 

In the diurnal tidal band, the results shown in Table 7 show that the satellite derived tide 
coefficients generally agree well with the oceanographic estimates from Schwiderski. The largest 
exception is the previously mentioned anomaly in S\ C 31 . The 0.4 cm correction in the K x C 21 
coefficient is the source of the approximately 10 mas term in the Lageos inclination residuals shown 
in Figure 2. In terms of in-phase and out-of-phase parts of the -S'] tide, the Lageos result agrees with 
Schwiderski for the out-of-phase part (C->i cos £ 21 ), but is 0.6 cm smaller for the in-phase part 
(-C 21 sin £ 21 ). The effect of changing the Earth’s free core nutation (FCN) period from 460 days to 
430 days as is indicated by analysis of nutation observations [Herring et al., 1991] should reduce the 
effective K x coefficient by 0.34 cm [Zhu et al., 1991], so the observed Lageos value is in better 
agreement with the smaller FCN period. Uncertainties in the actual ocean tide part of the K t signal 
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limits the usefulness of a more quantitative assessment of the FCN period deduced from the Lageos 
orbit analysis. The Starlette result from Cheng et al. [1990] is also more consistent with the lower 
FCN period. 

In the semi-diurnal tide results given in Table 8, two interesting points can be noted. First, both 
the Lageos and Starlette analyses confirm that the Schwiderski value for A/ 2 degree 2 order 2 
amplitude and phase requires a substantial correction. The out-of-phase part of this tide is 
C22 cos £22 and is the largest contributor to the tidal deceleration of the Moon’s mean motion 
[Cheng et al., 1992]. The contribution to the Moon’s h based upon the Lageos M 2 result is 
-20.35 arcsec century -2 , while that due to the Starlette result is -19.27. The Schwiderski value gives 
-15.15. Results in Marsh et al. [1990] derived from a combination of satellites agree well with the 
Lageos value, and the higher energy disspation obtained by all of the satellite results matches the 
observed secular deceleration results from lunar laser ranging [Dickey et al., 1990] much better than 
the Schwiderski value. Also, both the Lageos and Starlette results for S 2 C22 and e 22 include the 
contribution of the atmosphere as discussed in Chapman and Lindzen [1970]; however, the 
difference between these results and those of Schwiderski are substantially smaller than the predicted 
atmospheric effect. 

Although the adjustment of ocean tide parameters can remove many signals from the Lageos 
orbital element residuals, this parameterization can not handle the non-periodic variations of the low 
degree Stokes coefficients that are caused by non-tidal mass transfer in the atmosphere and oceans. 
Nor can they completely remove the effects of mismodeled nongravitational effects on Lageos. 

The analysis of the residuals in the ascending node is shown in Figure 3. The large curvature of 
the residuals using the nominal dynamical model in Figure 3a is caused by a combination of 
unmodeled secular change in / 2 an 18.6- year period variation due to uncertainties in the 
response on the inelastic Earth to long period tidal forces. The separation of these two effects, even 
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using more than 15 years of data, is still problematical due to extreme sensitivity of the results to 
unmodeled 7 2 variability at frequencies of less than 0.2 cycles per year [Eanes and Watkins, 1991]. 

Leaving the issue of long period and secular variability of J<i for further analysis, the node 
signals at frequencies greater that 0.2 cycles per year can be studied. Figure 3b shows the derivative 
of the node residuals of Figure 3a. Signals in this derivative, which will be referred to as the node 
excitation, are caused primarily by variations in J 2 and J 4 [Cheng et al., 1989]. Figure 3c shows the 
power spectrum of the time series of node excitations. The two largest peaks in the spectrum are 
labeled S a and S sa and represent annual and semiannual variability respectively. These anomalously 
large peaks are the source of the large adjustments of the degree 2 order 0 long period ocean tide 
coefficients discussed above. Although the semiannual peak is adequately removed by adjusting the 
S sa ocean tide, many attempts at removing all of the annual power by adjusting the S a tide fall short 
of this goal. The reason for this failure can be understood by the results, shown in Figures 3d and 3e, 
of a complex demodulation of the node excitation time series at the annual frequency followed by a 
bandpass filter. The results show that both the amplitude and the phase of the annual variation of J 2 
show substantial interannual variability. The amplitude shows variations of ±50 percent about the 
mean value of 100 mas year" 1 and the phase changes by ±30 degrees or about 1 month. The phase 
definition of Figure 3e is different from that of Table 3 by 93 degrees. A node excitation of 
100 mas year* 1 corresponds to an equivalent ocean tide amplitude of 1.98 cm or a J 2 variation of 
amplitude 24 x KT 11 . 

The stochastic nature of the observed node excitation precludes the possibility of solving for 
parameters in a deterministic model of the J 2 variability. Future studies of the Lageos results to the 
predictions using meteorological data must go beyond comparisons of the mean annual and 
semiannual terms to test the coherence of the two time series at all frequencies. Preliminary results 
indicate that substantial coherence exists in the frequency range from 1 to 4 cycles per year 
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[Gutierrez et al., 1991; Eanes and Watkins, 1991]. 

The most puzzling feature of the orbital element residuals of the Lageos long arc is the large 
signal in eccentricity (e) and argument of perigee (GO). These two elements are closely related and 
for a nearly circular orbit a change of variables to the nonsingular pair, e cos © and e sin co, simplifies 
the analysis [Yoder et al., 1983; Cheng et al., 1989]. The right-hand side of the differential equation 
for the complex quantity P = ecos©- i e sin© is referred to as the eccentricity vector excitation, 
HV Variability in the odd zonal harmonics cause variations with the same spectrum in the real part 
of 4V Errors in the odd degree diurnal and semi-diurnal ocean tide coefficients of order 1 and 2 
respectively cause variations in both the real and imaginary parts of 'F/>. 

Figures 4a and 4b show the real and imaginary parts of 'F/> using the nominal dynamical model 
of the long arc. The real part of *¥p is dominated by an annual variation while the imaginary part is 
dominated by a variation with a period of 560 days, period of the S ] tide perturbation and the 
interval of time required for one rotation of the Lageos orbital plane with respect to the Sun. These 
two features explain the anomalous estimates of the S a C\q and S\ C\\ tide coefficients in Tables 3 
and 4. But as mentioned above, the amplitudes of these terms are too large to be explained by tide 
effects alone, and Starlette analyses do not agree with these large values. This leads to the 
conclusion that this anomaly must be the result of a nongravitational acceleration of unknown origin 
acting on the Lageos orbit. The thermal imbalance and asymmetric albedo models (discussed later) 
that explain most of the observed Lageos drag-like acceleration and inclination slope do not seem to 
be adequate to remove the anomalous motion of the eccentricity vector. Until more is known about 
the source of the acceleration we must restrict ourselves to an empirical study which focuses on the 
form of the perturbations it causes. 

More insight into the nature of these perturbations is achieved by the complex demodulation of 
the real part of shown in Figure 4c and 4d. The amplitude of the annual term has a mean of 
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70masyear~ 1 corresponding to an equivalent ocean tide of 11 cm. From 1978 through 1987, this 
amplitude seems to decrease slightly and shows a modulation period of 3 years which is the beat 
period between the S\ and S a periods and which physically corresponds to the period of variation of 
the inclination of the Lageos node with respect to the ecliptic or the motion of the Lageos node with 
respect to the equinox (the K x period). In late 1988, a very large anomaly began which resulted in 
the amplitude tripling to more than 200 mas year' 1 during 1989. More recent analysis indicates that 
the amplitude also reached this same level in 1991. The spectrum of W P indicates that the 560 day 
term in the imaginary part of also exists in the real part, and together they are equivalent to a 
prograde oscillation in 'Fp with an amplitude of 40 mas year” 1 . Radial or transverse accelerations in 
near resonance with the orbital motion of Lageos are required to explain these observed trends. The 
required amplitude of the once per revolution acceleration has a peak value of 200 x 10” 12 ms’ 2 if it 
is in the transverse direction and twice this size if it is radial. This is about 1 percent as large as the 
direct solar radiation acceleration on Lageos. The eccentricity perturbations during the most 
anomalous year, 1989, reach 0.3 x 10 -6 , corresponding to error in the radial component of the Lageos 
position of more than 3 meters. Note that short arc fits to the Lageos data will accommodate this 
effect into the adjustment of the orbit initial conditions and lead to a much reduced signal in the 
range residuals. 

General Relativity 

One of the factors to be considered in the analysis of precise laser ranging data is the inclusion of 
general relativity in the data reduction procedure. The relativistic treatment of the near- Earth 
satellite orbit determination problem involves relativistic equations of motion, corrections to the 
measurement model, and time transformations. The problem can be formulated in a solar system 
barycentric reference frame or a geocentric reference frame. As a result of the generalized principle 
of equivalence, the main relativistic effects on a near-Earth satellite should be due only to the 
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Schwarzschild field of the Earth itself [Ashby and Bertotti, 1984]. Analysis of laser range data to 
Lageos in both reference frames has verified that the geocentric frame is adequate for Earth satellite 
applications [Ries et al M 1988; Huang et al., 1990]. The modeling for the barycentric frame is 
described in Martin et al. [1985], but Ries et al. [1988] used Lageos range data to demonstrate that 
an additional correction to the relativistic barycentric equations of motion is necessary to properly 
model the oblateness of the Earth’s gravitational field. 

The largest dynamical effect of general relativity on a satellite orbit is the well-known precession 
of perigee. For Lageos, the perigee precession is approximately 9 mas/day, an effect easily observed 
in the Lageos perigee residuals if not modeled. Relativity theory also predicts several small effects, 
including a change in the mean motion of the satellite and a precession of the longitude of the node 
due to the angular momentum of the rotating Earth, the Lense-Thirring effect [Lense and Thirring, 
1918]. In addition, there is a precession of the orbit plane due to the effect of geodesic (or de Sitter) 
precession [de Sitter, 1916]. This dynamical effect is not due to the Earth’s mass, but rather to the 
motion of the Earth through the Sun’s gravitational field. It is the result of the choice of a geocentric 
reference frame which is non-rotating with respect to the barycentric frame instead of a truly inertial 
geocentric reference frame, since the latter is difficult to realize in practice [Huang et al., 1990]. 
Observation of these small precessions by means of the analysis of the Lageos inclination and node 
residuals is presently not possible because of the uncertainties in the even zonal harmonics of the 
Earth’s gravitational field. It has been proposed, however, that the effects of the uncertainties in the 
even zonals can be eliminated by placing a Lageos-type satellite in an orbit identical to Lageos but 
with an inclination supplementary to Lageos [Ciufolini, 1986, 1989, Tapley et al., 1989]. 

The Earth’s Gravitational Coefficient 

The value of the gravitational coefficient of the Earth (GM) is an important parameter in the 
determination of the scale of the coordinate system realized by satellite observations [Zielinski, 
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1981]. The best values of GM are generally accepted as those determined by observing the influence 
of this parameter on the motion of near-Earth satellites, and the satellite best suited for this purpose 
is the Lageos satellite. Lageos was designed to minimize the effects of nongravitational forces, and 
the gravitational effects of all but the longest wavelength components of the Earth’s geopotential are 
greatly attenuated because of its high altitude. Thus the accuracy of the modeling of the forces on 
the Lageos satellite is more accurate than for any other satellite. In addition, the laser ranging 
measurements to Lageos are of very high accuracy. 

In 1985, a solution for GM using only laser ranging to Lageos produced a value of 
398600.434 ± 0.002 km 3 /sec 2 [Smith et al., 1985]. At the same time, a value for GM equal to 
398600.440 km 3 /sec 2 was determined from eight years of laser ranging to Lageos by UTCSR 
[Tapley et al., 1985], The uncertainty in the UTCSR estimate was reported subsequently as 
0.002 km 3 /sec 2 [Chovitz, 1987]. More recently, UTCSR reported a solution for GM obtained from a 
3-year fit to Lageos laser ranging, and also from a multi-satellite solution, to be 
398600.4405 + 0.001 km 3 /sec 2 [Ries et al., 1989]. The relativistic effects appropriate to the 
geocentric frame were modeled, where the coordinate time is equivalent to the current definition of 
Terrestrial Dynamical Time (TDT) [Huang et al., 1990]. 

In the initial determination, the Lageos laser range data were processed with a small but 
significant error in one of the range corrections. Optical tests on the Lageos-II satellite, which was 
built by the Agenzia Spatiale Italiana (ASI) to be a replica of the Lageos satellite, prompted a review 
of the tests conducted on Lageos. It was discovered that the value for the correction for the offset 
between the Lageos center-of-mass and the effective reflecting surface should be about 251 mm 
[Fitzmaurice et al., 1977] rather than the adopted 240 mm value [McCarthy, 1989]. In a new 
UTCSR determination of GM using the corrected center-of-mass offset, a value of 
398600.4415 km 3 /sec 2 in TDT units has been obtained, with an estimated uncertainty (l-o) of 
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0.0008 km 3 /sec 2 [Ries et al., 1992], 

While the direct effects of general relativity were taken into account in the data analyses 
described above, there is an indirect effect on the units being used that must be considered. The 
primary effect of general relativity on time is that coordinate time in different reference frames may 
run at different rates. Because the International Astronomical Union (IAU) definition of the time 
coordinate of the solar system barycentric reference system requires that only periodic differences 
exist between Barycentric Dynamical Time (TDB) and TDT, the spatial coordinates in the 
barycentric frame have effectively been rescaled [Misner, 1982; Hellings, 1986; Guinot and 
Seidelmann, 1988]. Thus the value of GM would be 398600.4418 km 3 /sec 2 in SI units and 
398600.4356 km 3 /sec 2 in TDB units. Recommendations have been made to modify the IAU 
definition for coordinate times to eliminate the rescaling, which would result in the units of length 
remaining SI units in all reference frames [Guinot, 1991]. 

Nongravitational Forces 

The dominant nongravitational force on Lageos is radiation pressure. The radiation may come 
directly from the Sun, it may be sunlight reflected by the Earth, or it may be infrared radiation that is 
emitted by the Earth. Temperature imbalances on the satellite itself can lead to ’photon thrusts’. 
Finally, there is expected to be some atmospheric drag even at the altitude of Lageos. Other forces, 
such as perturbations by the Earth’s magnetic field, interplanetary dust, the solar wind, or Poynting- 
Robertson drag, are not expected to have significant effects on the Lageos orbit [Rubincam, 1982; 
Ciufolini, 1989]. 

Solar and Earth Radiation Pressure 

Although the shape of the Lageos satellite is simple, modelling the effect of radiation pressure 
directly from the Sun is not trivial. A conical shadow model (umbra and penumbra) for the Earth 
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and Moon is generally used, but the numerical integration step-size is usually much larger than the 
duration of partial shadowing, and the effect is essentially a discontinuity in the solar radiation 
pressure force. Implementation of a procedure to more accurately account for the effect of 
shadowing is underway [Feulner, 1990]. While some orbit error may be attributable to the current 
method of integrating across the shadow boundary, experiments varying integrator step-size and 
shadow radius indicate that this is not the cause of the anomalous along-track accelerations observed 
in the Lageos orbit. 

It is unclear how much the modelling of the effects of Earth radiation pressure have benefitted 
the Lageos long arc fits. The effects are significant, and it has been possible to estimate the average 
reflectivity of the Earth from a Lageos long-arc, but it is likely that the orbital effects could be 
absorbed to some degree by other dynamical model parameters. The UTCSR model for Earth 
radiation pressure numerically integrates the heat and diffusely reflected sunlight from the visible 
disk of the Earth [Knocke et a 1., 1988]. The average and seasonal variations in the Earth’s albedo 
and infrared emission are included in a second-degree zonal representation. The temporal and 
geographic variations of the Earth’s albedo are expected to depart considerably from the zonal 
averages, so, at best, the Earth radiation pressure model only represents the long period effects. It is 
possible that some of the short period variations in the Lageos along-track acceleration are due to 
sunlight reflected diffusely and specularly by the Earth [Anselmo et al., 1983]. 

Surface Forces and the Lageos Spin Vector 

After subtracting most of the known forces acting on the Lageos satellite, there still remained a 
significant along-track acceleration which reduces the semi-major axis by approximately 1 mm per 
day. The mean values of the anomalous acceleration (or drag) determined every 15 days by UTCSR 
for the first 14.1 years of the previously discussed data set are plotted in Figure 5. The mean of the 
acceleration over the entire arc is about 3.5 picometer sec -2 with fluctuations that are sometimes as 
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large as the mean. The largest variations (spikes) occur when Lageos is experiencing eclipsing of 
the Sun by the Earth, although every eclipsing interval does not necessarily generate a spike. Thus 
both the mean and the variations required explanation. 

Earth Yarkovsky 

A thermal drag model, a variant of the Yarkovsky effect, has been proposed by Rubincam [1987, 
1988] which is able to account for much of the observed average along-track acceleration. The 
infrared radiation from the Earth is absorbed by the laser retroreflectors, and, because Lageos is 
spinning rapidly, the heat distribution is uniform longitudinally but not latitudinally. This creates a 
temperature imbalance between the Lageos northern and southern hemispheres, generating a thrust 
along the spin axis as the heat is re-radiated. The proposed thermal drag model accounts for about 
70% of the observed drag, and the remainder most likely consists of a combination of neutral particle 
drag and charged particle drag [Rubincam, 1990]. This model also predicts periodic variations about 
the mean with frequencies of once and twice per node revolution of the Lageos orbit (1050 and 525 
day periods). 

The effect of thermal drag due to Earth heating was simulated for the first 12.4 years of the 
Lageos mission. The 15-day averages of the accelerations generated by the nominal Rubincam 
[1988] model, augmented by 1 picometersec -2 to account for neutral and charged particle drag, 
compared well with the observed accelerations if the spikes were ignored. It was noted, however, 
that the modeled drag diverged from the observed drag in the last part of the 12.4 year arc in both 
amplitude and phase The deviation appeared to suggest that the spin axis of Lageos is evolving from 
its original orientation (where the agreement between the nominal model and observed drag is good) 
and becoming aligned with the Earth’s poles. A similar conclusion is reached by Rubincam [1990]. 
Rubincam [1987] suggested that the spin axis of the Lageos satellite should eventually align itself 
with the Earth’s spin axis. This has been confirmed by Bertotti and less [1991], who have analyzed 
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the effect of gravitational and magnetic torques on the satellite as the eddy current dissipation 
reduces the spin rate. In the UTCSR model, a model for the spin axis evolution was determined 
empirically so that the thermal drag fits the observed drag better. In Figure 6, it can be seen that the 
modified model is able to maintain good agreement in amplitude and phase of the periodic variation 
throughout the 12.4 year arc. A few direct measurements of the actual Lageos spin vector are clearly 
needed, since this effect, and the effects described below, depend strongly on the orientation of the 
spin axis. 

Examination of the orbit inclination residuals provide additional evidence supporting the thermal 
drag model. In the UTCSR Lageos long arcs where the thermal drag was not modeled, the residual 
errors in the orbit inclination have exhibited a slope of 1.3 to 1.5 mas year -1 . The slope in the 
Lageos inclination residuals had been a concern, since there seemed to be no reasonable force which 
could generate this particular signal. The effect of a co-rotating atmosphere at Lageos altitude was 
considered, but it was found that even if one assumed that 100% of the drag was atmospheric drag, a 
co-rotating atmosphere could account for less than 10% of the inclination slope. By incorporating 
the ’Earth Yarkovsky’ model in the latest long arc, the inclination slope is only 0.3 mas year -1 . It is 
convincing to note that the thermal drag model is able to explain much of the average drag, the 
variation at the 1050 and 525 day periods, and most of the peculiar inclination slope. If the 
magnitude of the thermal drag force was increased to about 90% of the observed drag, then the 
remaining inclination slope would be explained also. Any other explanation for the inclination slope 
is still unknown. 

Solar Yarkovsky 

A similar thermal thrust mechanism, due to heating by the sun, has been proposed to account for 
at least part of the large spikes that occur only during eclipse seasons [Rubincam, 1982; Slabinski, 
1988; Afonso et al., 1989; Scharroo et al., 1991; Ries et al., 1991]. Since the spin axis orientation is 
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essentially fixed with respect to the Sun during an orbit, the hemisphere of the Lageos satellite which 
is experiencing summer (tilted towards the sun) will be warmer than the opposite hemisphere at 
every point in the orbit. Ignoring the Earth’s heat and the effect of shadowing, this temperature 
difference will be essentially constant, and there is no significant net effect on the orbit. During 
shadowing, however, the solar thermal force does not average out, and there is a net along-track 
acceleration. 

A model for this effect was included in the force model for Lageos, and the estimated drag was 
compared to the observed drag. The maximum magnitude of the ’Solar Yarkovsky’ acceleration 
(80 picometersec -2 ) was empirically chosen to give peaks with amplitudes comparable to the 
observed drag. Afonso et al. [1989] and Scharroo et al. [1991] independently analyzed the expected 
surface temperatures and arrived at accelerations with similar magnitudes. The spin axis orientation 
in the UTCSR model was based on the modified spin axis model when calculating the ’season’ for 
each hemisphere. Unfortunately, it has not been possible to choose a set of parameters for this 
model which could predict the correct series of peaks. 

Asymmetric Reflectivity 

An additional mechanism is proposed that, when combined with the Solar Yarkovsky effect, 
appears to account for the spikes observed during the eclipse seasons. If it is assumed that the two 
halves of the satellite do not have the same effective reflectivity, then there will be an asymmetry in 
the solar radiation pressure on the satellite that, averaged over the spin period, will be along the 
direction of the spin vector. Like the solar Yarkovsky effect, there is no significant orbital effect 
except during eclipse seasons. The cause of the reflective asymmetry is not known. Scharroo et al. 
[1991] speculates that the two halves of the satellite may have been finished differently and finds that 
the northern hemisphere of Lageos need only be l/70th more specular than the southern one. Ries et 
al. [1991] suggests that the non-symmetrical distribution of the infrared comer-cubes may be the 
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cause. There are three infrared cubes in the southern hemisphere and only one in the northern 
hemisphere, and since they appear opaque at optical wavelengths, the northern hemisphere is likely 
to be more reflective than the southern hemisphere. Whatever the source of the asymmetry, 
augmenting the thermal forces with an imbalance in the reflectivity leads to a model that is capable 
of generating the accelerations similar to those observed on the Lageos satellite over the first 12.4 
years as shown in Figure 7. The agreement, however, is still not perfect, which illustrates that the 
models are still deficient in some respects. 

Lageos Spin Vector 

As the Lageos spin vector evolves with time, the models described above may become less 
reliable. Analysis of the along- track accelerations estimated in the latest Lageos long-arc, which 
incorporated the UTCSR surface force models described above, indicates that the agreement after the 
first twelve years is degrading. It is critical for the modeling of these surface forces that 
measurements of the spin axis orientation are obtained. Eventually, however, the spin rate will slow 
enough to result in chaotic behavior. When this occurs, it is not clear whether the thermal forces will 
average out and thus be attenuated, or become more significant and more difficult to model. 

DEFINITION OF THE TERRESTRIAL REFERENCE FRAME 

Using the short arc adjustments as discussed earlier enables not only an improved fit to the laser 
ranges as measured by range residual rms, but also allows a frequency domain separation of orbit 
error from kinematic effects in the residuals, namely those due to the positions of the ground 
tracking sites. These sites, being tied to the surface of the Earth, are forced to have a diurnal 
variation in the inertial frame. In the satellite frame, the variation differs from once per revolution 
by one cycle per sidereal day a signal which is easily detectable and separable from other effects. 
The most significant remaining orbit error at this frequency is due to error in the order 1 geopotential 
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coefficients, but the effect is considerably smaller than that of the kinematic signal of the site 
position errors at the few millimeter level when averaged over many revolutions. 

The rigid rotations of the tracking network, that is, that part of the site motions common to all the 
sites in rotation, are, by definition, polar motion and Earth rotation, collectively referred to as Earth 
orientation parameters (EOP). Clearly, if both the complete tracking network and all EOP’s are 
adjusted simultaneously, a singularity results, and the arbitrary orientation of the network aliases into 
the EOP’s. This problem is solved primarily by adjusting the network positions once over a 
relatively long span of time while adjusting the EOP’s frequently, such as once per day. Under these 
restrictions, only the mean values of the EOP series are inseparable from the orientation of the 
network. This last singularity is generally removed by the application of one of several types of 
constraints, including the fixing of specific fiducial site(s), or the fixing of EOP’s on one reference 
day [Smith et al., 1991; Robertson et ah, 1990]. An alternative method is to apply an analytic 
constraint equation that forces the "net rotation" of the tracking network adjustments with respect to 
a nominal set to be zero, and they are thus absorbed by the EOP series [Bender and Goad, 1979]. 
Such a constraint allows continuity of the combined tracking network and EOP series, referred to as 
the terrestrial reference frame (TRF), over time as more data from additional sites is acquired. The 
implementation of this method is under study at UTCSR, although for the present, an ad hoc 
realization of this approach is used by constraining through a priori covariance both the EOP and 
tracking sites. This covariance is chosen to be sufficiently loose to allow short period variation in 
the EOP’s without allowing long-term variations or drift. For the results presented in this paper, the 
a priori uncertainty on each site coordinate was 1 meter, 10 mas on x and y pole position, and 0.7 ms 
on UT1. A priori correlations were set to zero. By not fixing any single site, the terrestrial reference 
frame is freed from the vagaries, in either data quality or quantity, of a particular site. 
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The estimation of UT1 using SLR, or any satellite technique, is complicated by the fact that long 
period modeling of the motion of the longitude of the ascending node of the orbit(s) is at present 
impossible due to the stochastic variation in the mass distribution of the atmosphere, and to 
unmodelled surface forces. The relative size of these perturbations depend upon tha satellite design 
and orbit characteristics. For Lageos, the variability in zonal mass distribution is the more 
significant. Thus, UTCSR currently produces UT1 estimates which are constrained to a VLBI 
solution at long periods, but which are increasing independent as the frequency increases. This 
technique was first described in Robertson et al. [1983], using a Gaussian filter with full width at half 
maximum of 90 days, and yielded rms agreement of 0.7 ms. By contrast, using an improved 
Vondrak smoothing technique tying to JPL SPACE90 [Gross et ah, 1991], and greatly improved 
observations and analyses from both techniques, the rms agreement in 1990 was less than 0.07 ms 
[IERS Ann. Rep., 1990]. 

Inclusion of Site Velocities 

When the span of observations becomes long enough, site velocities may need to be estimated, 
particularly for sites in deforming regions where the a priori plate model may not be accurate. The 
inclusion in the adjusted parameter set of site velocities adds an additional singularity to the above 
discussion, namely that between the mean slope in the EOP’s and the net rotation rate derived from 
the combined velocities of the sites in the TRF. This is resolved in an analogous manner to that 
described above, through fixing fiducial velocities, two reference EOP days, or explicit constraint 
equations. The current implementation at UTCSR involves constraint through a priori covariance, 
and through the adjustment only of those sites with considerable velocity departures from the 
nominal Minster and Jordan AMO-1 velocity [Minster and Jordan, 1978]. Alternatively, many sites 
can be adjusted, and the resulting EOP series can be detrended with respect to an a priori series for 
continuity. This approach is used when many sites velocities are adjusted for studies of global 
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tectonics. 


Site Position, Velocity, and Earth Orientation Solutions from UTCSR 

With the above introduction as background, the several types of solutions for geodetic 
parameters produced at UTCSR can be discussed. The first type of solution is performed annually 
and represents an updated version of the UTCSR TRF. This solution spans the entire Lageos 
mission from 1976 through the end of the reporting year, and adjusts mean site coordinates, 
velocities for selected sites, and a time series of EOP’s, most recently using 3 day resolution. Since 
the series is reported to the International Earth Rotation Service and forms a significant portion of 
the combined IERS EOP series and International Terrestrial Reference Frame (ITRF), it is important 
that these series be expressed in a well defined reference frame, and hence not all sites, even those 
with sufficient data, have adjusted velocities. In the most recent such solution, SSC(CSR)91L03, 
seven sites in tectonically deforming regions had velocity adjustments: 7110 Monument Peak 
(USA), 7109 Quincy (USA), 7838 Simosato (Japan), 7939 Matera (Italy), 7907 Arequipa (Peru), 
7097 Easter Island (Chile), and 7123 Huahine (French Polynesia). 

The EOP series associated with this solution are analyzed for the combined IERS series by the 
IERS Central Bureau at the Paris Observatory. For the 1990 report, the EOP(CSR)91L03 series was 
assessed at roughly the 0.6 mas level in x and y polar motion, and 0.07 ms in UT1, for the period 
1986-1990 [IERS Annual Report, 1991]. This was among the most accurate of any series reported 
that year. A typical polhode of recent polar motion as determined from UTCSR Lageos and IRIS 
VLBI is presented in Figure 8. The considerable progress made by all space geodetic techniques is 
evident when the history of such intercomparisons is recalled. Robertson et al. [1983] presents 5 
mas agreement in pole position between UTCSR SLR and IRIS VLBI. Robertson et al. [1985] 
presents agreement during 1984 — 1985 of 2 mas in pole position. Comparisons in BIH and IERS 
annual reports from 1986 to 1990 trace the accuracies to their current submilliarcsecond level. 
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Interestingly, Carter et al. [1984] also presented one of the first comparisons of EOP derived from 
space geodetic techniques and that predicted from measurements of atmospheric angular momentum 
(AAM), and demonstrated clearly the significant El Nino event of 1983 and associated peak in 
excess length of day. The seasonal discrepancy between space geodetic measurements and AAM 
was also noted, and it was correctly conjectured to be largely due to the effects of the atmosphere 
(stratosphere) above 50 mb. 

A second type of solution produced on an annual basis adjusts all sites and velocities (with 
sufficient data) for the basis of assessing global tectonic models. Roughly 30-40 sites have strong 
velocity solutions [Watkins et al., 1989; Smith et al., 1991]. Relative motions for sites in Europe 
derived from such solutions is presented in Figure 9. Alternatively, to test internal consistency of the 
site solutions, independent solutions of few month duration can be determined for each site, and the 
velocities inferred from these time series [Watkins, 1990]. Such a solution for the baseline length 
between Monument Peak, CA, on the Pacific plate, and Quincy, CA, on the North American plate is 
presented in Figure 10. The slope of the best fit line is the well known —25 mm/yr adjustment with 
regard to the RM2 predicted velocity, and the rms scatter about the best fit line of 7 mm is 
remarkably small for such a long observation period. 

The accuracy of the epoch site coordinates can best be measured through comparison with an 
independent technique with collocated observations. A careful comparison of SLR solutions 
computed at UTCSR and VLBI solutions computed at GSFC for sites occupied primarily in 1988 
resulted in rms agreements of 15, 21, and 22 mm in x, y, and z, respectively [Ray et al., 1990]. These 
numbers include the uncertainties in each solution and the uncertainties in the survey ties connecting 
the observing monuments of each technique. Analysis of the chi-square per degree of freedom 
indicated that each techniques formal uncertainties should be scaled by approximately 2, resulting in 
uncertainties of less than 10 mm in each component for the best observed sites. An extension of the 


1.30 



Ray et al. paper analyzing the results of site positions determined using significantly more data, and 
requiring mapping with estimated velocities for both techniques is presented by Himwich et al. [this 
issue]. 

By computing successive seven parameter transformations between the several month solutions 
and the nominal TRF, the translational origin or geocenter position can be determined. The most 
significant proposed source of motion between the solid Earth, on which the tracking sites reside, 
and the center of mass of the solid Earth-atmosphere-hydrosphere is mass redistribution in the ocean. 
Figure 1 1 demonstrate the geocenter history in 15 day intervals since 1987. The rms about the mean 
is 12 mm in X, 9 mm in E, and 25 mm in Z. The source of the long period trend in X is unclear, but 
it is statistically significant (rms after removal of best fit line is 9 mm), and may represent true long 
period geocenter motion. Full analyses of the geocenter time series is the subject of a forthcoming 
manuscript. 

A third type of solution is performed each week, and is referred to as the UTCSR Earth 
Orientation Rapid Service. This service involves the gathering of recent Lageos data, and the 
computation of the EOP’s in near real time, and provides prompt reports to the United States Naval 
Observatory and the IERS for their rapid service Bulletins A and B. This solution has been 
performed and reported weekly since 1982. Site coordinates in this operational service are held fixed 
to the values from a previous annual solution and not readjusted, hence the EOP’s from the 
operational service are perhaps a few tenths of a milliarcsecond less accurate than those of the annual 
solutions. 

Future Improvements to SLR Derived Geodetic Parameters 

A number of improvements lie on the horizon for the determination of geodetic parameters using 
the SLR technique. The first of these is improved temporal resolution of EOP’s obtained through 
Kalman filtering of the orbit parameters. Because the SLR data are obtained each day, although 
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some days may be sparse, information on the orbit and geodetic parameters is available on a daily 
basis, although some knowledge of the spectral power of the unmodeled orbit excitations and high 
frequency EOP excitation is required for optimal performance. The improved temporal resolution of 
the space geodetic EOP series, which may also be possible with the use of GPS, will aid in the 
understanding of momentum exchange between the components of the Earth system, particularly 
between the atmosphere and solid Earth, at short periods. Although a manuscript on the application 
of Kalman filtering of the Lageos orbit and the resulting high frequency EOP series is in preparation, 
preliminary results are presented in Figure 12, which demonstrates the UTCSR Lageos 1 day x and y 
pole positions plotted with GSFC processing of the high density, high quality VLBI ERDE 
(Extended Research and Development Experiment) period in the fall of 1989. Both series are 
differenced from a smoothed Lageos nominal series. The overall rms agreement in each coordinate 
is 0.7 mas, although slightly better during the heart of the campaign (October). In addition, 
determination of the diurnal and semidiurnal variations in EOP’s induced primarily by ocean tides 
has been demonstrated at the 20 microarcsecond level or better through analysis of a recent span of 
high quality Lageos data from 1987-1991 [Watkins et aL, 1991]. 

Another exciting development is the use of additional geodetic satellites such as Lageos-2 and 
the Soviet Etalon-1 and Etalon-2. Lageos-2, an identical twin of the current Lageos- 1, will be 
launched in fall 1992 into an orbit similar to that of Lageos- 1 but with an inclination of 52 degrees. 
Assuming adequate tracking is available for both satellites, the use of observations to both satellites 
in a simultaneous solution can reduce the time on site for SLR systems by approximately a factor of 
two while retaining the same accuracy as currently obtained. The Etalon spacecraft, launched into 
orbits of the Glonass spacecraft, similar to those of the Global Positioning System, can also provide 
reduced time on site, but can also assist in the extension of the period over which the SLR derived 
UT1 estimates can remain independent from those of VLB! This is made possible by the 
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attenuation at the Etalon altitude of the stochastic variations in even zonal harmonics of the Earth’s 
gravity field which drive the satellite node. The response of the nodal longitudes of the Etalon 
spacecraft is roughly 10 times less to a given excitation in even zonal harmonics than those of 
Lageos-1 or Lageos-2. Thus the current limit of approximately 50 days over which SLR can obtain 
independent measurements of UT1 may be extended considerably. 

CONCLUSIONS 

Satellite laser ranging has matured over the last decade into one of the essential space geodesy 
techniques. It has demonstrated centimeter site positioning and millimeter per year velocity 
determinations in a frame tied dynamically to the mass center of the solid Earth-hydrosphere- 
atmosphere system. Such a coordinate system is a requirement for studying long term eustatic sea 
level rise and other global change phenomena. Earth orientation parameters determined with the 
coordinate system have been produced in near real time operationally since 1983, at a relatively 
modest cost. The SLR ranging to Lageos has also provided a rich spectrum of results based upon the 
analysis of Lageos orbital dynamics. These include significant improvements in the knowledge of 
the mean and variable components of the Earth’s gravity field and the Earth’s gravitational 
parameter, The ability to measure the time variations of the Earth’s gravity field has opened as 
exciting area of study in relating global processes, including meteorologically derived mass transport 
through changes in the satellite dynamics. Finally, new confirmation of General Relativity has been 
obtained using the Lageos SLR data. 

With the launching of Lageos-2 and Lageos-3, the future decade holds significant promise for 
substantially improving the accuracy and applicability of the SLR technique. 
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Fig. I. Evolution of the laser ranging network tracking Lageos from 1976 to 1991. 

Fig. 2. Analysis of inclination residuals from Lageos Long Arc 9107. (a) Residuals using the 
nominal dynamical model, (b) Residuals after adjusting dynamical model parameters, (c) First 
derivative of the residuals in (a), (d) Power spectrum of the time series in (c). 

Fig. 3. Analysis of longitude of the ascending node residuals from Lageos Long Arc 9107. (a) 
Residuals using the nominal dynamical model, (b) First derivative of the residuals in (a), (c) Power 
spectrum of the time series in (b). (d) Amplitude and (e) phase of the annual variation obtained by 
complex demodulation of the time series in (b). 

Fig. 4. Analysis of the eccentricity vector excitation residuals from Lagesos Long Arc 9107. The (a) 
real and (b) imaginary parts of the eccentricity vector excitation 0F/>). The (c) amplitude and (d) 
phase of the annual variation of the real part of 'Vp obtained by complex demodulation of the time 
series in (a). 

Fig. 5. Observed along-track acceleration for Lageos. 

Fig. 6. Modified Earth Yarkovsky model compared to observed along-track acceleration during the 
first 12.4 years of the Lageos mission. 

Fig. 7. Combined models for Earth and Solar Yarkovsky, asymmetric reflectivity and atmospheric 
drag compared to observed accelerations. 

Fig. 8. Polar motion obtained from UTCSR Lageos and IRIS VLBI, 1987-89. 

Fig. 9. Baseline rates 

Fig. 10. Variability of the baseline length from Monument Peak to Quincy, 1982-90. 

Fig. 11. Estimates of the translations along the x, y, and z axes of the instantaneous origin of the 
terrestrial reference system with respect to the geocenter, 1987-91. 

Fig. 12. Kalman smoothed polar motion at a 1-day resolution from Lageos compared to GSFC VLBI 
results during the 1989 ERDE Campaign. 
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Figure 5. Observed along-track acceleration for Lageos 



Figure 6. Modified Earth Yarkovsky model compared to observed along-track 
acceleration during first 12.4 years 
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Figure 7. Combined models for Earth and Solar Yarkovsky, asymmetric reflectivity 
and atmospheric drag compared to observed accelerations 
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Figure 9. Relative baseline rates for selected European sites 













Figure 10. Baseline length variability. Monument Peak to Quincy, 1982-1990 
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Figure 12. Lageos 1 day Kalman filtered polar motion compared to GSFC VLBI 
during the 1989 ERDE Campaign 




